#include<iostream>
using namespace std;
#include<cstdlib>
#include<cmath>
#include<fstream>

/*
 * random(...)	产生在 [a,b] 区间内，以 func(x) 为概率密度的随机数，在使用这个函数之前，可以用srand()函数对随机数种子进行初始化
 * double a
 * double b			随机数所在区间 [a,b]
 * double (*func)(double x)	随机数的相对概率密度
 * double fmax			相对概率密度的上界
 */
double gnr_rand(double a, double b, double (*func)(double x), double fmax){

	double x = (double)rand() / RAND_MAX * (b-a) + a; // 在 [a,b] 之间均匀分布
	double f = (double)rand() / RAND_MAX * fmax; // 在 [0, fmax] 之间均匀分布
	while( f > func(x) ){ // x 有 func(x) 的相对概率被接受
		x = (double)rand() / RAND_MAX * (b-a) + a;
		f = (double)rand() / RAND_MAX * fmax;
	}
	return x;
}

double sample(int n, double a, double b, double (*func)(double x), double fmax){

	double y=0;
	for(int i=0;i<n;i++){
		y += gnr_rand(a, b, func, fmax);
	}
	return y/n;
}

/*
 * count_random(...) 统计随机数rand出现在 [a,a+h), [a+h, a+2h), ... 中的频率，记录在count中
 */
void count_random(double a, double h, int * count, double rand){

	count[ (int)((rand-a)/h) ] ++;
}

double f(double x){//[-3:3]之间均匀分布
	if(x >= -3 && x <= 3) return 1.0/6;
	else return 0;
}

double NormDistr(double mu, double sigma, double x){
	return 1 / sqrt(2*M_PI) / sigma * exp( - (x-mu)*(x-mu)/2/sigma/sigma );
}

int main(){

	int i,N=1000000;
	double a=-3, b=3, fmax = 1.0/6;
	double h=0.01, x;
	int n = (b-a)/h;
	int *count = new int [n];
	for(i=0;i<n;i++) count[i] = 0;
	for(i=0;i<N;i++){
		x = sample(30, a, b, f, fmax);
		count_random(a, h, count, x);
	}
	ofstream fp("test_LLNandCLT.txt");
	for(i=0;i<n;i++){
		fp<< a+ (i+0.5)*h <<"\t"<<(double)count[i]/N/h<<"\t"<< NormDistr(0, sqrt(3.0/30), a+ (i+0.5)*h )<<endl;
	}
	fp.close();
	delete [] count;

	return 0;
}
